/*
 * Compare if two matrices are identical, i.e. abs(a[i][j] - b[i][j]) < EPSILON.
 * It can be used to test only the upper triangle.
 */
#include <stdio.h>
#include <unistd.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>

#define EPSILON 0.001

char delim[] = "\t \r\n";

int main(int argc, char **argv) {
   FILE *f1, *f2;
   char *line, *word;
   size_t bufsize;
   ssize_t len;
   int ncol,nrow;
   int i,j;
   float x1,x2;
   int upper = 0;

   if (argc < 3) {
      fprintf(stderr,"Usage: %s <matfile1> <matfile2> [upper]\n",argv[0]);
      exit(1);
   } 

   if (argc > 3) {
      if (!strcmp(argv[3],"upper")) {
         upper = 1;
      }
   }

   if ((f1 = fopen(argv[1],"r")) == NULL) {
      fprintf(stderr,"Cannot open file %s\n",argv[1]); 
      exit(1);
   }

   if ((f2 = fopen(argv[2],"r")) == NULL) {
      fprintf(stderr,"Cannot open file %s\n",argv[2]); 
      exit(1);
   }

   
   // determine the number of columns
   line = NULL;
   if ((len = getline(&line,&bufsize,f1)) < 0) {
      fprintf(stderr,"Cannot read a line in %s\n",argv[1]); 
      exit(1);
   }
   ncol = 0;
   while ((word = (char *) strtok(line,delim)) != NULL) {
      ncol++;
      line = NULL;
   }
   printf("The number of columns = %d\n",ncol);
   rewind(f1);
   free(line);

   i = 0; j = 0;
   while (!feof(f1)) {
      if (fscanf(f1,"%f",&x1) == EOF) 
         break ;
      fscanf(f2,"%f",&x2);
      //if (x1 != x2) {
      if ( !upper || ((upper) && (i < j)) ) {
         if (fabs(x1 - x2) > EPSILON) {
            printf("(%d,%d): %f <> %f\n",i,j,x1,x2);
         }
      }
      j++; 
      if (j == ncol) {
         i++;
         j = 0;
      }
   }
   nrow = i;
   if ((nrow != ncol) || ((j % ncol) != 0)) {
      fprintf(stderr,"Not a square matrix\n");
      exit(1);
   }
   printf("%d lines read\n",nrow); 
   fclose(f1);
   fclose(f2);
}
